****************************************************************************
**		Authors: 	Or Tuttnauer & Liran Harsgor
**		Purpose: 	Produce Figure 5 in Kedar, Harsgor & Tuttnauer (JOP)  
**		input:		KHT_countrylevel.dta
*****************************************************************************

log using "figure 5.log", replace

set more off
use "KHT_countrylevel.dta", clear

gen ctr = substr(cyear, 1, 3)
keep country cyear ctr enps enpv lnmeddm sr smd mmm fused vsr vlmeddm

*** A. GENERATE DIFFERENCES ***

** generate simulated difference **
* rerun our model
reg enps enpv lnmeddm sr vlmeddm vsr smd mmm fused, robust
predict yhat
save original_data.dta, replace 

* add uncertainty *
scalar rss=e(rss)  			/*sum of squared residuals*/
scalar dfr=e(df_r)			/*degrees of freedom*/
scalar sig=sqrt(rss/dfr)	/*sigma square hat*/
matrix msig=(sig)
drawnorm eps, n(54) means(0) sds(msig)
gen ystar=yhat+eps

* run baseline model on simulated DV *
reg ystar enpv lnmeddm    vlmeddm     smd mmm fused, robust
	matrix b=e(b)
	scalar b1=b[1,1]  		/* first coefficient (enpv) */
	scalar b3=b[1,3]		/* interaction coefficient */
gen pred1=b1+b3*lnmeddm

* run our model on simulated DV *
reg ystar enpv lnmeddm sr vlmeddm vsr smd mmm fused, robust
	matrix b=e(b)
	scalar b1=b[1,1]  		/* first coefficient (enpv) */
	scalar b4=b[1,4]		/* first interaction coefficient */
	scalar b5=b[1,5]		/* second interaction coefficient */
gen pred2=b1+b4*lnmeddm+b5*sr

gen diff=pred1-pred2
save bootstrap.dta, replace

** B. NOW PRODUCE 4999 ADDITIONAL DIFFS FOR EACH COUNTRY 
 
local i=2
forvalues i = 1(1)4999 {
	use original_data.dta, clear
		drawnorm eps, n(54) means(0) sds(msig)
	gen ystar=yhat+eps

	* baseline model
	reg ystar enpv lnmeddm    vlmeddm     smd mmm fused, robust
		matrix b=e(b)
		scalar b1=b[1,1]  		/* The first coefficient (enpv) */
		scalar b3=b[1,3]		/* The interaction's coefficient */
	gen pred1=b1+b3*lnmeddm

	* our model
	reg ystar enpv lnmeddm sr vlmeddm vsr smd mmm fused, robust
		matrix b=e(b)
		scalar b1=b[1,1]  		/* The first coefficient (enpv) */
		scalar b4=b[1,4]		/* The first interaction's coefficient */
		scalar b5=b[1,5]		/* The second interaction's coefficient */
	gen pred2=b1+b4*lnmeddm+b5*sr

	gen diff=pred1-pred2

	save temp.dta, replace
	use bootstrap.dta, clear
	append using temp.dta
	save bootstrap.dta, replace
}


bysort cyear: egen meandiff=mean(diff)

** generate CIs ** 
bysort cyear: egen p25=pctile(diff), p(2.5)
bysort cyear: egen p975=pctile(diff), p(97.5)

egen country_unique = tag(cyear)

*** Generate graph ***
graph twoway rcap p25 p975 sr if country_u==1 & smd==0 & mmm==0, legend(off) scheme(s1mono) ///
yline(0) ytitle("{&Delta} permissibility") xtitle("Seat ratio") ylabel(,labsize(vsmall)) ///
|| scatter meandiff sr if country_u==1 & smd==0 & mmm==0, mcolor(black) msymbol(D) ///
msize(vsmall) || scatter p975 sr if country_u==1 & smd==0 & mmm==0, msymbol(i) ///
mlabel(ctr) mlabsize(vsmall) mlabposition(12)
graph save "figure 5.gph", replace 

log close
